#########################################################################################
#########################################################################################
####Non-state Atrocities in Capital Cities - Mechanism 3 (Diplomatic Activity Plots) ####
#########################################################################################
#########################################################################################
library(foreign)
library(plyr)
library(DataCombine)
library(doBy)
library(countrycode)
library(ggplot2)

# Set working library
setwd("~/Data/Mechanism 3/")
# Read in data on insurgent atrocities in the capital (1996-2009) data
cap.list <- read.csv("caplist.csv")
# Read in the data on total diplomatic activity for 1995
dip.dat.1995 <- read.dta("Dip_ex_95.dta")

####################################################
###Create plotting data from raw activity dataset###
####################################################
#Change the "other" category to NAs so that index focuses on 1=chargé d'affaires; 2=minister; 3=ambassador; 0=no evidence of diplomatic exchange
dip.dat.1995$DR_at_1[dip.dat.1995$DR_at_1==9]  <- NA
#Average the level of diplomatic activity in a given country's capital in 1995 by collapsing the dydadic data by country
dip.sum.1995 <- summaryBy(DR_at_1 ~ ccode1, FUN=c("mean"), data=dip.dat.1995, na.rm=TRUE)
names(dip.sum.1995) <- c("ccode", "DR1995")
#Merge with the total number of insurgent atrocities in each country's capital (total, 1996-2009)
j <- join(dip.sum.1995, cap.list)
#Remove country's with missing information
j <- na.omit(j)
#Create a log version of the total number of insurgent atrocities in the capital
j$logcapat <- log(j$capat+1)

#######################
###Correlation Plots###
#######################
###Full data
#Correlation
cor(j$DR1995, j$logcapat)
#Trendline
lm.pred <- lm(logcapat ~ DR1995, data=j)
##Confidence intervals
#Use predicted values and extract var-covar matrix
pp.pred <- predict(lm.pred)
V <- vcov(lm.pred)
X <- model.matrix(~ DR1995, data=j)
#Calulate stanard erors
se.fit <- sqrt(diag(X %*% V %*% t(X)))
#Combine all to a dataframe and create standard errors
predframe <- with(j, data.frame(DR1995,
                                logcapat=pp.pred,lwr=pp.pred-1.96*se.fit,upr=pp.pred+1.96*se.fit))

#Create baseline graph
p <- ggplot(j, aes(x = DR1995, y = logcapat))
#Plot full graph
pdf("empbat.pdf")
p +  labs(x="Average Level of Diplomatic Relations, 1995", y="Capital Cell Years Experiencing Insurgent Atrocities, 1996-2009 (Log)") +
  geom_point()+ 
  #scale_y_discrete(limits=c(0, 3.5)) +
  annotate("text", x = 2, y = 3, label = "Correlation = 0.12") +
  geom_line(data=predframe, size=0.5, color="red") +
  geom_ribbon(data=predframe,aes(ymin=lwr,ymax=upr),alpha=0.3)
dev.off()


###Atrocity outliers removed
#Subset data to only cities that experienced less than 10 insurgent atrocities within the 1996-2009 period
k <- subset(j, capat<10)
#Correlation
cor(k$DR1995, k$logcapat)
#Trendline
lm.pred <- lm(logcapat ~ DR1995, data=k)
##Confidence intervals
#Use predicted values and extract var-covar matrix
pp.pred <- predict(lm.pred)
V <- vcov(lm.pred)
X <- model.matrix(~ DR1995, data=k)
#Calulate stanard erors
se.fit <- sqrt(diag(X %*% V %*% t(X)))
#Combine all to a dataframe and create standard errors
predframe <- with(k, data.frame(DR1995,
                                logcapat=pp.pred,lwr=pp.pred-1.96*se.fit,upr=pp.pred+1.96*se.fit))

#Create baseline graph
p <- ggplot(k, aes(x = DR1995, y = logcapat))
#Plot full graph
pdf("empbatno.pdf")
p +  labs(x="Average Level of Diplomatic Relations, 1995", y="Capital Cell Years Experiencing Insurgent Atrocities, 1996-2009 (Log) -- Outliers Removed") +
  geom_point()+ 
  #scale_y_discrete(limits=c(0, 10)) +
  annotate("text", x = 2, y = 2.25, label = "Correlation = 0.17") +
  geom_line(data=predframe, size=0.5, color="red") +
  geom_ribbon(data=predframe,aes(ymin=lwr,ymax=upr),alpha=0.3)
dev.off()

